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The Canonical Transformation theory of Yanai and Chan [J. Chem. Phys. 124, 194106 (2006)] 
provides a rigorously size-extensive description of dynamical correlation in multireference problems. 
Here we describe a new formulation of the theory based on the extended normal ordering procedure 
of Mukherjee and Kutzelnigg [J. Chem. Phys. 107, 432 (1997)]. On studies of the water, nitrogen, 
and iron-oxide potential energy curves, the Linearised Canonical Transformation Singles and Dou- 
bles theory is competitive in accuracy with some of the best multireference methods, such as the 
Multireference Averaged Coupled Pair Functional, while computational timings (in the case of the 
iron-oxide molecule) are two-three orders of magnitude faster and comparable to those of Complete 
Active Space Second-Order Perturbation Theory. The results presented here are greatly improved 
both in accuracy and in cost over our earlier study as the result of a new numerical algorithm for 
solving the amplitude equations. 



I. INTRODUCTION 



lated contributions as we describe below. 
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While nondynamic correlation between electrons es- 
tablishes the qualitative features of chemical bonding, it 
is the accurate description of dynamic correlation, asso- 
ciated with the short-range cusp behaviour of the wave- 
function, which is necessary to obtain quantitative agree- 
ment with experiment. Starting from a suitable refer- 
ence function, the exponential ansatz provides an accu- 
rate and economical description of dynamic correlation. 
For example, in systems that are qualitatively described 
by a single determinant reference, Coupled Cluster (CC) 
theory paired with a large basis set yields predictions 
with chemical accuracy [H, H, [H • However, for the many 
chemical problems which require a multireference char- 
acterisation, a practical theory for dynamic correlation 
with the desirable qualities of the exponential ansatz - 
size-extensivity, chemical- accuracy, and moderate com- 
putational cost - has yet to be widely established. 

In an earlier article [3] we presented a Canonical Trans- 
formation (CT) theory which is based on an exponential 
ansatz, is rigorously size-extensive, and which may easily 
be combined with any multireference starting wavefunc- 
tion. In the form implemented in that work the compu- 
tational cost is 0(a re ), where a is the number of active 
orbitals and e is the number of external orbitals. In cal- 
culations of bond-breaking potential energy curves, the 
Linearised Canonical Transformation Doubles (L-CTD) 
theory performed significantly better than multireference 
perturbation theory, and obtained the accuracy of Cou- 
pled Cluster Single Doubles (CCSD) at the equilibrium 
geometry across the entire potential energy curve. Our 
work was directly motivated by the Canonical Diagonal- 
isation theory of White @ although there are earlier re- 
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The purpose of the current work is to improve on our 
initial contribution in several areas. A central feature 
of the Canonical Transformation theory is the use of an 
operator decomposition, both to close the infinite expan- 
sions associated with an exponential ansatz and to reduce 
the complexity of the energy and amplitude equations 
that arise when working with a complicated reference 
function. In our earlier work, we introduced a cumulant- 
type operator decomposition by analogy to the cumulant 
decomposition of density matrices found in reduced den- 
sity matrix theories @, 0, H, H, [HI [H[. However, this 
choice of operator decomposition is not unique and here 
we explore an alternative operator decomposition, with 
some formal advantages, that is based on the concept 
of extended normal orde ring as introduced by Mukherjee 
and Kutzelnigg [T2, [H, Il4j . Indeed, examination of the 
articles by these authors shows that they anticipated the 
utility of their results in multireference correlation theo- 
ries, and in this context our current theory is in part a 
realisation along such directions. 



A second focus of this work is to investigate in detail 
the behaviour of the Canonical Transformation theory in 
a variety of chemical problems. For example we study, 
with a range of basis sets, the bond-breaking potential en- 
ergy curves of water, nitrogen, and iron-oxide and com- 
pare our results against state-of-the-art multi-reference 
configuration interaction and perturbation theories. In 
addition, we examine numerically the size-extensivity 
and density-scaling properties of the Canonical Trans- 
formation energies. The results in the present study are 
much improved over our earlier work, in large part be- 
cause of improvements we have made to our numerical 
algorithms, and we describe in detail the numerical as- 
pects of efficiently implementing and converging the CT 
equations. 
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II. CANONICAL TRANSFORMATION THEORY 
A. Recapitulation 

In multireference problems we divide the orbitals into 
active orbitals, which describe the nondynamic correla- 
tion and external orbitals which describe the dynamic 
correlation. The external orbitals may further be di- 
vided into core and virtual orbitals; core orbitals are 
those which remain doubly occupied in all the reference 
configurations. 

We will assume that a reference wavefunction is 
available that describes the nondynamic correlation in 
the problem. This may be obtained, for example, from a 
Complete Active Space Self- Consistent Field (CASSCF) 
calculation that exactly correlates electrons within the 
active orbitals [HI, [l6| . Alternatively, and especially for 
larger active spaces, a Density Matrix Renormalization 
Group wavefunction may be used [13, HH • We then in- 
corporate the remaining dynamic correlation on top of 
the reference wavefunction via an exponential oper- 
ator that generates excitations between the active and 
external spaces yielding 

* = e A * (1) 

We will be concerned with a unitary formulation, where 
A> = —A. The excitations are understood to be both of 
external and semi-internal form 

A = At(a? al) + A$(a$ - <%,) + KjKj -<&) + - 

(2) 

where ijk. . . denote active indices, abc. . ., external in- 
dices, a° = a\ai, afj = a\a b ajai, and the summation 
convention is assumed. For example, the first two terms 
are the usual external single and double excitations, while 
the third term (with three active indices) is a semi- 
internal single excitation, which captures the coupling 
between singles relaxation in the active space and singles 
excitation to the external space. 

In a related picture, we can also view e A as generating 
an effective canonically transformed Hamiltonian H that 
acts only in the active space, but which has dynamic 
correlation folded in from the external space, where 

H = e- A He A (3) 
= (4) 

The exponential ansatz combined with a multirefer- 
ence wavefunction as shown in eqn. (TT|) has a long 
history and we necessarily can only give an incomplete 
account here. Such an ansatz is used in some forms of 
multi-reference coupled cluster theory (MRCC) as dis- 
cussed in the review by Paldus and Li [19(. In particular, 
an early example of a complete theoretical scheme for a 
related multi-reference coupled cluster method was given 
by Mukherjee in Ref. [12fl . While CC theory is usu- 
ally formulated in terms of similarity rather than canon- 
ical (i.e. unitary) transforms, unitary exponentials have 



previously been explored in a multi-reference setting by 
Freed et al [2(|, Kirtman et al [21], and Simons et al 
[22| . We mention also the single-reference unitary cou- 
pled cluster work by Kutzelnigg [H, [24[ , Bartlett et al 
[H H| H3, and Pal [H, HI]. The general concept of ef- 
fective Hamiltonians and canonical transformations is of 
course very old, dating back to van Vleck [30]. We note 
in particular some modern theories that emphasize an 
effective Hamiltonian language similar to our own such 
as the Effective Valence Hamiltonian theory of Freed I20l| 
and the Generalized van Vleck theory of Kirtman [2l| . 
As recognised by Freed, the folding in of dynamic cor- 
relation into the active-space effective Hamiltonian is a 
form of renormalisation transformation. This picture was 
pursued by White in his theory of Canonical Diagonalisa- 
tion [5j, and as described previously, this is the primary 
precursor to our work. 

In the exponential ansatz of single-reference coupled 
cluster theory, the commutativity of the excitation oper- 
ators in the single-reference form of A, i.e. A — Afaf + 
Afjafj, allows the Baker-Campbell-Hausdorff expansion 
of H to terminate at low-order for low-particle rank in A. 
The difficulty in working with the multireference expo- 
nential ansatz arises from the non-commuting excitations 
in the multi-reference form of A in eqn. @ , which leads 
to a non-terminating expansion for the effective Hamil- 
tonian H . (In fact this difficulty already arises if we use 
a unitary e A with the single-reference form of A). 

In our earlier Canonical Transformation (CT) theory 
we introduced a new route to a tractable and com- 
putationally efficient formulation for the multireference 
ansatz (TTJ). Starting from the Baker-Campbell-Hausdorff 
expansion of the exact effective Hamiltonian, 

H = H+[H,A] + ^[[H,A],A] +... (5) 

we replace each commutator by an approximate decom- 
posed commutator, to yield an approximate effective 
Hamiltonian 

#1,2,... =H+[H,A] lt2i ... + ^[[H,A] 1A ...,A]] 1A ... + ... 

(6) 

Each subscript denotes a decomposition, and the num- 
bers 1,2... denote the particle ranks of the operators 
that remain after the decomposition. Note that if all par- 
ticle ranks were included in the decomposition (i.e. the 
subscripts ranged from 1, 2, ... n, where n is the number 
of particles), then eqns. (O and © would be identi- 
cal. If in addition to including all particle ranks in eqn. 
^ A contained up to n-body excitations, then the CT 
ansatz (JXJ) would be exact in the sense of full configura- 
tion interaction, and indeed eqn. ^ would hold exactly. 
The two relevant approximations thus arise from restrict- 
ing the excitations in A (wavefunction ansatz) as well as 
the form of the operator decomposition (operator ansatz) 

M 

As an example, let us consider the linearised CT single 
and doubles theory (L-CTSD) introduced in our earlier 
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work. Here A is restricted to contain only one- and two- 
particle excitations as in eqn. (J2J), and we restrict all 
decomposed commutators to contain at most one and 
two-body operators (i.e. subscripts 1,2). Since [H, A] 
generates a three-body operator, this requires some de- 
composition of a three-body operator into lower body 
operators. We proposed an explicit decomposition into 
one- and two-body operators based on an analogy to the 
cumulant decomposition of density matrices, 

aZ => 9(7? A aZ) 12(<tf A 7 « A <) (7) 

where in the above A denotes an antisymmctrisation over 
all upper and lower indices with an associated factor of 
l/(p) 2 , i.e. 1/36 in the above case, where p is the particle 
rank of the original operator. Here we will present the 
explicit steps leading to the above decomposition. Our 
notation follows closely that of Kutzelnigg and Mukherjee 
[3ll ] . Recall that the cumulant decomposition provides a 
way to rewrite reduced density matrices 7 in terms of 
products of cumulants A, via 

7? = K) = K (8) 
7 P J = (a™) = \ p J+^- 1 ?~/ 9 s (9) 

iTu = «> = KZ + iX - itKu + iX 
-iX + iX - iX + iXu - iXu + iX 
+i p s ihu - iWtil + liilil - 
Wsihl - i r s ih p u (10) 

For the three-particle density matrix, by dropping the 
three-particle cumulant A^, and substituting the ex- 
pressions © and ^ in (flTj)) . we obtain an approximate 
decomposition in terms of one- and two-particle density 
matrices only 

pqr p qr p qr , p qr 

I stu 7s ^itu it 7stz ' luist 

q pr , q pr q pr . r pr r pr _■_ r pr 
7s7iu 1 it 7su 7?z7si ' ^snu it 1 su ' 7u7si 

-2(7?« - 7?7[7* + iHYu ~ iHlu 

= 9(7rA7 t l r )-12(7?A 7t 9 A 7 D (11) 

To obtain our operator decomposition, we simply re- 
placed expectation values in the above terms by the 
corresponding operators, i.e. j^t ~ * a ^t an d 7? - * of, 
yielding eqn. ([7]). Note that by construction, the expec- 
tation value of the operator decomposition reproduces 
the three-particle density matrix cumulant decomposi- 
tion (fTT|) . 

By using this decomposition recursively i.e. by con- 
structing the double commutator by first using the de- 
composed single commutator [H, A]\ y 2 as in eqn. ©, the 
full effective Hamiltonian Hi 2 at the L-CTSD level con- 
tains only one and two-body operators. Evaluation of the 
energy then only requires the one- and two-particle den- 
sity matrices of the reference function. As discussed in 
our initial work, this fulfils one of the criteria for an ef- 
ficient multireference theory, namely, we do not need to 



explicitly manipulate the complicated reference function. 
From a different perspective, the canonical tranforma- 
tions can also be viewed as providing a parametrisation 
of a two-particle density matrix theory. Recently, such 
connections have been explored from a different direction 
by Mazziotti [12, [H[ and while interesting, we shall not 
dwell further on these matters here. 

We call the above formulation a linearised theory, be- 
cause the operator decomposition is applied at the first 
commutator. Then, at the L-CTSD level the energies 
and amplitudes are evaluated via 

E= (¥„ \&i,2 |*o) (12) 

= (*o|[ J ff 1 , 2) <-<] 1 , 2 |*o) (13) 

= <M#i ia ,a#-a?Ji, 2 |tt o > (14) 

= <*o|[ffi,a,<- 0^1,21*0) (15) 

The resulting computational cost of the theory is 
0(a 2 e ), and is thus comparable to that of a single- 
reference coupled cluster calculation. 

B. Accuracy of the operator decomposition 

As presented above, the accuracy of the Canonical 
Transformation theory rests on the accuracy of opera- 
tor decomposition, given at the L-CTSD level by eqn. 
([?))■ However, although our operator decomposition was 
chosen so that its expectation value would reproduce the 
density matrix cumulant decomposition, this choice is 
not unique. For example, we could add to the r.h.s. of 
eqn. ([7]) any term with vanishing expectation value with 
^0 and still preserve the correspondence with the den- 
sity matrix cumulant decomposition (jlip . This simply 
reflects the fact that a decomposition for expectation val- 
ues (i.e. the cumulant decomposition) does not contain 
sufficient information to specify a corresponding operator 
decomposition. 

In our earlier work, we examined the accuracy of the 
operator decomposition through a perturbative analysis 
of CT theory starting from a single determinantal wave- 
function tyf) and using a single- reference single-doubles 
excitation operator A = Af(af - a} a ) + A$(a$ - a l j b ). 
This analysis showed that the L-CTSD theory was ac- 
curate through third-order in the fluctuation potential 
W = H — F where F is the Fock operator i.e. 

(y D \H\V D ) 
=<* 'D\Si t2 \$ D ) + 0(W 4 ) 

=(* D \H + [H, A]i, a + [[H, A} li2 ,A} li2 \^ D ) + 0(W A ) 

(16) 

However, consider what happens if we use the more 
general multireference form of A in eqn. ((2|) that in- 
cludes semi- internal excitations such as Af*(afj° — a% k ), 
together with a single reference wavefunction |^£>). Such 
excitations should not contribute as they destroy the 
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single reference wavefunction, and thus all expecta- 
tion values of exact commutators containing only semi- 
internal excitations, e.g. (V D \[H,A$(a$ - a% k )]\V D ), 

(* D | [[H, A# (a% - , A b ^ n « - a%)] \ * D ) must van- 
ish. However, using the cumulant-based operator decom- 
position |7|) we find that although the expectation value 
of the first commutator {* D \[H,A$(a$ - <4)]i, 2 |*i?), 
correctly vanishes, it does not do so for the second com- 
mutator. Non-vanishing terms arise e.g. from 

(*n\[[H,A$(a% - a« )]i, 2 , - <4«)]l^> (17) 

Writing H and the two A operators as g^g^gg, o'o'ov, 
v'o'oo respectively, using g,o,v to denote general, occu- 
pied, and virtual indices respectively, we can see a non- 
zero contribution arising from 

i — i i — i 

(^oKsVss) (oVov) (uVoo)!*!)) ^ (18) 

where the underbracket denotes contraction and the over- 
bracket denotes a replacement by a density matrix in the 
operator decomposition. 

In a multireference situation, we use the same ex- 
tended excitation operator A (with semi-internal excita- 
tions) through the entire potential energy surface, even 
when the underlying reference wavefunction is largely of 
a single reference nature, as is sometimes the case near 
the equilibrium geometry. Thus the above deficiency of 
the cumulant operator decomposition for single reference 
wavefunctions motivates us to examine other possible de- 
compositions, as we describe now. 



III. EXTENDED NORMAL ORDERING 



A. Normal ordering for a Multireference 
wavefunction 



Normal ordering provides a standard way to decom- 
pose an operator into a sum of zero-, one-, two- and 
higher body contributions that are ordered with respect 
to a given vacuum. In many-body theory it is common 
to use normal ordering not with respect to the physical 
vacuum, but rather with respect to a single determinant 
state or Fermi vacuum. With respect to the Fermi vac- 
uum, normal ordering of the operators a p s , a pq s , a pq ^ yields 



~aP s =aP s - S p n s 



(19) 
(20) 



°? s t ~ 5 p s n s a q t - 8 q t n t a p s 
+ S P n t a q s + S q n s a P + S Pq n p n q 

~pqr _ pqr _ ™ qr cq pt , cr qp 
a stu — a stu °s n P a tu T °s n q a ru + °s n r a tu 

+ S p n p af u - 5 q n q a p ; u + 6 r t n r a pq u + 5 p u n p aZ 

+ 5q u n q a pr t - 5 r u n r a pq + n p n q S p Ja r u + n p n r 5 pr u a q 

+ n q n r SZa p - n p n q S Pq a r t - n p n q S pq a r s - n p n r 5 p s r t a q u 

- n p n r S P ^a q - n q n r bf s a p u - n q n r 8 qr u a p t - n p n q n r 8 p s Z 

(21) 



where the tilde represents operators normal-ordered 
w.r.t. the Fermi vacuum (quasi-particle operators), n p 
is the occupation number (0 or 1) of the p-th orbital, 
and S% = SPSj - 5 p 5«, S P Z = + S P S q Jl + S>5«% - 

8 p t 8 q s 8 r u - 5P6q t 6 r s - 5 p 8 q J r t . Note that all normal-ordered 
operators (other than the "zero-body" constant term) 
yield a vanishing expectation value with the Fermi vac- 
uum, e.g. (af) = 0. If we arc interested in a state which 
is well approximated by the Fermi vacuum, the higher- 
particle rank quasi-particle operators such as a pq ^ are less 
relevant to its properties than the lower-rank ones, since 
they represent multiple simultaneous excitations away 
from the state. Thus the Fermi-vacuum normal ordering 
presents a natural way to approximate high-particle rank 
operators in terms of simpler lower-body terms by simply 
neglecting the high particle-rank quasi-particle operators 
that appear in the normal ordered form. For example, to 
approximate a s q u in terms of one- and two-body opera- 
tors alone, we would neglect a pq ^ in cqn. (f2"Tj) . 

In the Canonical Transformation theory, however, we 
are often interested in reference states which cannot 
be represented well by any Fermi vacuum. Recently, 
Mukherjee and Kutzelnigg proposed an elegant gener- 
alisation of normal-ordering w.r.t. such multireference 
states [lH, [l4| . By examining the form of the above 
normal-ordering equations when rotated into an arbitary 
one-particle basis, they arrived at the generalised rela- 
tions 



■r. 



(22) 

Kt + 1>! + i q t~< - ft~< - + il? 

a?? + 4( 7 ? A~a q ) + 1 pq (23) 



1st 
s,1ZP V 



a stu * Is a tu Is a tu Is a tu It a su * It a su 
_ ~r~pq _ p~qr _ q~pr T ~pq pq~ r , pr ~q 
It a su Iu a ts Iu a st > Iu a st + lst a u T lsu a t 
, qr~p nn~r pq~r pr~a pr~a qr~ p 

+ 7ttX - lsu a t ~ lui a s - 1st a l - ltu a l - Its < 

= ~<7u + 9(7? A ~aZ) + 9( 7 ^ A <) + 7^ (24) 



Let us examine the physical meaning of the above ex- 
pressions, taking eqn. (|23p as an example. Here, we see 



that the original two body operator a pq is written in 
terms of an average over the reference state (the zero- 
body operator 7? t 9 ), a product of a one-body average 
with a one-body quasi-particle operator (the terms like 
7f5,j), and a two-body quasi-particle operator a p s q . The 
quasi-particle operators describe fluctuations about the 
reference, because just as in the usual form of normal or- 
dering, their expectation values with the reference vanish 
e.g. (a? ) = 0, (a p J) = 0. 



B. Application to Canonical Transformation 
Theory 

The extended normal ordering provides a systematic 
operator decomposition which is well suited to Canonical 
Transformation theory. At the linearised CTSD level, we 
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wish to decompose the three-body operators, arising from 
the commutator [H, A], into lower-body terms. We can 
do so by neglecting the effects of the simultaneous three- 
body fluctuations described by the operator aPJ^. For 
consistency, we should also remove the fully connected 
three-body cumulant Af^. First let us rewrite af^ in 
terms of a p s ,a p s 1 be rearranging eqn. (|24p , and substitut- 
ing in the cumulant decomposition of 7^ (I10p , we find 

„pqr ~pqr _ -,J>r fl «r _ qi r _ r\ , 
u stu — a stu Is l u tu It \ a u lu) ' • ' • 
ori 1 pa / r r , 

+ AX+7?A^ + ... 7 ?7^ + --- 
- 5SS + 9(7? A <) - 36( 7 f A 7? A <) + 9( 7 ^ A <) 
+24( 7 p A 7 « A <) - 9(tJ A 7 £) + A^ (25) 

Now dropping S,^ and Ag' u we obtain the extended 
normal-ordered decomposition, which we name the MK 
decomposition after Mukherjee and Kutzelnigg, 

<Z =*»(t? A <) - 36( 7 f A 7? A <) + 9( 7s ^ A <) 
+24( 7 PA 7t 9 A 7 ;;)-9( 7 ^A<) (26) 

Comparing the MK decomposition to our earlier 
cumulant-type decomposition ([7]) we see that they yield 
the same expectation value with the reference function 
*f> and thus differ only by terms whose expectation val- 
ues vanish. In addition to some different factors, the MK 
decomposition include additional operators: a constant 
term, and the term 7 ^ t 9 /\a r u . Computationally, both these 
terms are easily implemented without affecting the scal- 
ing of the original L-CTSD algorithm. 

To better understand the differences between the MK 
and cumulant-type (CU) decompositions, it is instructive 
to compare the two for a simpler example, namely, the 
decomposition of the two-particle operator a p s q . These 
are 

a^2( 7 J A a?) 

= i (7 ^ + 7t X- 7t X- 7 |a?) CU 

(27) 

a^ 7 ^- 7 <<) + 7t 9 K- 7 f) 

- lt(a q s ~ li) ~ l q sK - 7? ) + Ysll ~ 7?T2 MK 

(28) 

Here we see that the MK decomposition is expressed in 
terms of fluctuations e.g. a q s — 7 | in the presence of the 
field 7 f , while the cumulant decomposition involves the 
bare operators a q s directly. The neglected term a v s 1 in 
the MK decomposition has the conceptual meaning of a 
simultaneous two-particle fluctuation operator, and we 
consider this to be conceptually appealing. 

Returning to the earlier example that motivated our 
examination of alternative operator decompositions, let 
us now look at the normal-product decomposition of com- 
mutators involving semi-internal excitation operators, as 
in eqn. (|17[) . Starting from a single determinantal refer- 
ence, the extended normal ordering reduces to the usual 



TABLE I: Total energies of FCI and differences of various 
methods from FCI for the simultaneous bond breaking of H2O 
molecule with CAS(6e, 5o) and cc-pVDZ basis sets. The units 
are E h . The bond angle is fixed at <HOH = 109.57°. R e 
= 0.9 929 A. r s = 10~ 2 and r d = 10" 2 (described in Sec. 
I VII B|) were used in the L-CT calculations. See Ref. Q for 
the previous L-CT results. 





1 u 




■irie 


4K e 


FCI 


-76.23885 


-75.94558 


-75.91003 


-75.90872 


RHF 


0.21718 


0.37002 


0.57365 


0.67159 


CASSCF 


0.16299 


0.13196 


0.12302 


0.12259 


CAsF 1 2 


0.01330 


0.00843 


0.00848 


0.00852 


CASPT3 


0.00377 


0.00383 


0.00174 


0.00158 


MR-CI 


0.00556 


0.00378 


0.00296 


0.00290 


MR-CI+Q 


-0.00056 


-0.00053 


-0.00066 


-0.00068 


MR-ACPF 


0.00093 


0.00054 


0.00020 


0.00017 


MR-AQCC 


0.00231 


0.00150 


0.00102 


0.00098 


CCSD 


0.00384 


0.02248 


0.00967 


0.00200 


CCSDT 


0.00051 


-0.00238 


-0.04106 


-0.04973 


L-CTSD(CU) 


0.00029 


-0.00097 


-0.00171 


-0.00172 


L-CTSD(MK) 


-0.00077 


-0.00128 


-0.00192 


-0.00192 


previous L-CTD(CU) 


0.00219 


-0.00056 


0.00297 


0.00251 


previous L-CTSD(CU) 


0.00061 


-0.00358 


0.00301 


0.00287 



normal ordering with respect to a Fermi vacuum de- 
scribed by eqns (fl"$|) - (f2Tj) . Then, the operator decomposi- 
tion corresponds to dropping the three-particle normal- 
ordered operators aP^ in eqn. (f2"Tj) . By construction, 
the remaining normal-ordered operators e.g. af, q all have 
vanishing expectation value with the Fermi vacuum, and 
consequently using the MK decomposition, the expecta- 
tion values of all commutators of the form of eqn. (fT7|) 
with single determinant references vanish as they should, 
in contrast to the cumulant-type decomposition. 

Thus we see that the extended normal-ordered MK 
decomposition offers some conceptual and formal advan- 
tages over our earlier cumulant-type CU decomposition. 
Encouraged by these aspects, we have implemented this 
decomposition and we now proceed to the numerical re- 
sults. 



IV. CALCULATIONS 

A. Water and nitrogen potential energy curves 

We performed prototype multireference CT calcula- 
tions for the simultaneous bond breaking curve of the 
water molecule and the bond breaking curve of the 
nitrogen molecule. We chose these molecules to al- 
low a direct comparison with the results in our previ- 
ous paper with the CU decomposition. Here we 
have used a wider range of basis sets, including the 
cc-pVDZ and cc-pVTZ basis sets for water and 6-31G, 
cc-pVDZ, and cc-pVTZ basis sets for nitrogen [34l l35j. 
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TABLE II: Total energies of MR-CI+Q and differences of 
various methods from MR-CI+Q for the simultaneous bond 
breaking of H2O molecule with CAS(6e, 5o) and cc-pVTZ ba- 
sis sets. The units are E^. The bond angle is fixed at <HOH 
= 109.57°. R e = 0.9 929 A. t s = 10" 1 and r d = 10~ 2 (de- 
scribed in Sec. IVII Bfl were used in the L-CT calculations. 







2R e 


3-R e 


4R e 


MR-CI+Q 


-76.32847 


-76.01591 


-75.97484 


-75.97345 


RHF 


0.27679 


0.41697 


0.60895 


0.70500 


CASSCF 


0.22228 


0.18279 


0.16880 


0.16815 


CASPT2 


0.01545 


0.00915 


0.00904 


0.00911 


CASPT3 


0.00574 


0.00646 


0.00307 


0.00286 


MR-CI 


0.01005 


0.00761 


0.00611 


0.00603 


MR-ACPF 


0.00232 


0.00176 


0.00139 


0.00137 


MR-AQCC 


0.00467 


0.00353 


0.00280 


0.00277 


CCSD 


0.00742 


0.02995 


0.02724 


0.01999 


CCSDT 


-0.00055 


-0.00147 


-0.03965 


-0.04866 


L-CTSD(CU) 


0.00214 


0.00058 


-0.00081 


-0.00084 


L-CTSD(MK) 


0.00186 


-0.00004 


-0.00102 


-0.00103 



For assessment, we carried out calculations with state-of- 
the-art internally contracted multireference methods - 
second- and third-order perturbation theory (CASPT2 
and CASPT3) 1361 , 371 l38l . |39| |. configuration interaction 
(MR-CI) pfl l4ll. l42l| . the a posteriori size-extensivity cor- 



TABLE III: Total energies of FCI and differences of various 
methods from FCI for the bond breaking of N2 molecule with 
CAS(6e, 60) and 6-31G basis sets. The Units are Eh- t s = 
10 _1 and r d = 10 -2 (described in Sec. IVII Bfl were used in 
the L-CT calculations. See Ref. [|[ for the previous L-CT 
results. 





lA 


2A 


3A 


FCI 


-109.04667 


-108.85968 


-108.83905 


RHF 


0.21143 


0.55008 


0.85649 


CASSCF 


0.08551 


0.08623 


0.07472 


CASPT2 


0.01372 


0.00834 


0.00830 


CASPT3 


0.00558 


0.00769 


0.00409 


MR PT 


n nn9R8 

u.uuzuo 


u.uuouo 


n nn9i n 


MR-CI+Q 


-0.00012 


-0.00014 


-0.00016 


MR-ACPF 


0.00092 


0.00071 


0.00027 


MR-AQCC 


0.00133 


0.00125 


0.00069 


CCSD 


0.00685 


-0.00731 




CCSDT 


0.00122 


-0.05220 




L-CTSD(CU) 


0.00142 


-0.00165 


-0.00173 


L-CTSD(MK) 


0.00082 


-0.00187 


-0.00250 


previous L-CTD(CU) 


0.00510 


0.00447 




previous L-CTSD(CU) 


0.00646 


0.00112 





rected configuration interaction due to Davidson (MR- 
CI+Q) (Htm]], averaged coupled pair functional (MR- 
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FIG. 1: Energy differences E-E(FCI) for the simultaneous 
bond breaking of H2O molecule with CAS(6e,5o) and cc- 
pVDZ basis sets. 



FIG. 2: Energy differences E-E(MR-CI+Q) for the simulta- 
neous bond breaking of H2O molecule with CAS(6e,5o) and 
cc-pVTZ basis sets. 
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TABLE IV: Total energies of MR-CI+Q and differences of 
various methods from MR-CI+Q for the bond breaking of N2 
molecule with CAS(6e, 60) and cc-pVDZ basis sets. The units 
are E h . r a = 10" 1 and r d - 10~ 2 (described in Sec. IVII Bj) 
were used in the L-CT calculations. 





lA 


2A 


3A 


MR-CI+Q 


-109.22891 


-108.98376 


-108.96035 


RHF 


0.29907 


0.65317 


0.96627 


CASSCF 


0.18453 


0.19413 


0.18316 


CASPT2 


0.02243 


0.01558 


0.01616 


CASPT3 


0.00700 


0.00781 


0.00375 


MR-CI 


0.00926 


0.01161 


0.01011 


MR-ACPF 


0.00262 


0.00245 


0.00173 


MR-AQCC 


0.00419 


0.00465 


0.00374 


CCSD 


0.01112 


0.07424 




CCSDT 


0.00177 


-0.04382 




L-CTSD(CU) 


0.00118 


0.00024 


-0.00045 


L-CTSD(MK) 


0.00117 


0.00162 


0.00026 



ACPF) [H, Hfl and averaged quadratic coupled-cluster 
theory (MR-AQCC) [13] (both a priori size-extensivity 
modifications of configuration interaction), as well as 
single-reference coupled cluster calculations at the CCSD 
and CCSDT level (H, El])- Full configuration interac- 
tion (FCI) energies were also used for comparison where 
available. The CAS space for the multireference cal- 



TABLE V: Total energies of MR-CI+Q and differences of var- 
ious methods from MR-CI+Q for the bond breaking of N2 
molecule with CAS(6e, 60) and cc-pVTZ basis sets. The units 
are E h . t s = 10" 1 and r d = 10" 2 (described in Sec. IVII Bp 
were used in the L-CT calculations. 





lA 


2A 


3A 


MR-CI+Q 


-109.33774 


-109.05871 


-109.03045 


RHF 


0.36972 


0.70119 


1.00458 


CASSCF 


0.25475 


0.25035 


0.23571 


MR-CI 


0.01452 


0.01725 


0.01516 


MR-ACPF 


0.00363 


0.00324 


0.00234 


MR-AQCC 


0.00625 


0.00669 


0.00548 


CASPT2 


0.02399 


0.01094 


0.01195 


CASPT3 


0.00753 


0.01077 


0.00396 


CCSD 


0.01532 


0.09593 




CCSDT 


0.00021 


-0.03276 




L-CTSD(CU) 


0.00453° 


-0.00006 


-0.00051 


L-CTSD(MK) 


0.00249 


0.00162 


0.00035 



a) t s = 5 x 10 1 and Td — 2 x 10 2 were used because of 
convergence problems. 



culations was six active electrons in five active orbitals 
[denoted (6e,5o)] for the water calculations and (6e,6o) 
for the nitrogen calculations. The Is orbitals in O and 
N atoms were held frozen in all calculations. For the 
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FIG. 3: Energy differences E-E(FCI) for the bond breaking 
of N2 molecule with CAS(6e, 60) and 6-31G basis sets. 
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FIG. 4: Energy differences E-E(MR-CI+Q) for the bond 
breaking of N2 molecule with CAS(6e, 60) and cc-pVDZ basis 
sets. 
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L-CTSD calculations, we employed both the cumulant 
(CU) and normal-ordering (MK) operator decomposi- 
tions described in section [TTJ The internally-contracted 
multireference calculations were executed using molpro 
[11] , the CC calculations using tce [49[ in utchem [5(| , 
and the CT calculations using our own computer pro- 
gram. 

Tables [H HH HTH HVl and |V] present the errors in the 
total energies of various methods as measured from FCI 
or (in the larger basis sets) MR-CI+Q at several points 
across the potential curve. These errors are plotted in 
Figures [Ullllll and [U 

Comparing all the different methods, in the calcula- 
tions where FCI energies were available, MR-CI+Q pro- 
vided the smallest maximum absolute error (MAE) and 
non-parallelity error (NPE) and for this reason was used 
as the benchmark method when FCI energies could not 
be obtained. The general order of accuracy in terms of 
MAE from most to least accurate was MR-CI+Q « MR- 
ACPF « L-CTSD(CU), L-CTSD(MK) w MR-AQCC > 
CASPT3 « MR-CI > CASPT2. While the MAE of L- 
CTSD(CU) and L-CTSD (MK) was comparable to that 
of MR-ACPF and MR-AQCC, the NPE was larger; in 
the intermediate region the shapes of the curves some- 
what resembled the CASPT3 curve. In the equilibrium 
region, the L-CTSD energies were similar in accuracy to 
CCSDT. 

The MAE and NPE for the two CT operator decompo- 
sitions L-CTSD(CU) and L-CTSD(MK) are compared in 
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FIG. 5: Energy differences E-E(MR-CI+Q) for the bond 
breaking of N2 molecule with CAS(6e, 60) and cc-pVTZ basis 
sets. 



TABLE VI: Maximum absolute error (MAE) and non- 
parallelity error (NPE) of L-CTSD (CU) and L-CTSD (MK). 
The units are mE^. 

L-CTSD(CU) L-CTSD(MK) 





MAE 


NPE 


MAE 


NPE 


H 2 0/cc-pVDZ 


1.72 


2.01 


1.92 


1.11 


H 2 0/cc-pVTZ 


4.25 


5.10 


2.79 


3.82 


N2/6-3IG 


2.07 


4.08 


2.54 


4.74 


N 2 /cc-pVDZ 


3.35 


3.83 


5.09 


4.86 


Na/cc-pVTZ 


3.62 


4.46 


5.49 


5.35 



Table IVI1 We find that the two operator decompositions 
performed quite similarly in these systems, although the 
MAE of L-CTSD(CU) was slightly smaller. For compar- 
ison, we have also included the L-CTSD(CU) energies 
from our calculations in our earlier work [4J. We note 
that our new L-CTSD (CU) energies are significantly im- 
proved, particularly in the intermediate dissociation re- 
gion. This is a result of the new numerical algorithm, 
described in IVII B[ which allowed us to significantly re- 
duce the truncation of the operator manifold that we used 
in our previous work. However, the curves of the new L- 
CTSD in the figures are not completely smooth due to 
some remaining operator truncation effects in the numer- 
ical solution and removal of this non-smooth behaviour 
will be addressed in future work. 



Tablc lVril shows the spectroscopic constants of N2 com- 
puted by fitting the potential curves. Compared to the 
available FCI results in the 6-31G basis, MR-CI+Q once 
again came closest for all spectroscopic parameters (R e , 
w e , D e ) while the related MR-ACPF and MR-AQCC 
methods behaved very similarly to MR-CI+Q. Compar- 
ing CT against the other methods, different trends were 
observed for different quantities. For the dissociation en- 
ergies, we found that MR-CI+Q > MR-ACPF « MR- 
AQCC w L-CTSD(CU) > L-CTSD(MK) « CCSDT > 
CASPT3 > CASPT2 > CCSD. For frequencies, in the 
cc-pVDZ and cc-pVTZ basis L-CTSD was comparable in 
accuracy to MR-ACPF /MR-AQCC (though with errors 
in the opposite direction) and better than than those 
of CCSDT, while the equilibrium bond distances were 
less accurate than MR-ACPF /MR-AQCC though still 
comparable to CCSDT. L-CTSD (CU) and L-CTSD (MK) 
generally performed similarly, although the spectroscopic 
constants for L-CTSD (CU) with cc-pVTZ could not be 
obtained because of convergence problems at the fitting 
geometries. The small non-smoothness in the potential 
energy curves resulting from the numerical approxima- 
tions in solving the CT equations may also be a factor in 
the less systematic errors of the CT methods for R e and 
to e . 

Thus to summarise, the overall performance of L- 
CTSD(CU) and L-CTSD(MK) for these potential en- 
ergy curves was competitive with the best multirefer- 
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TABLE VII: Spectroscopic constants for N2 molecule by var- 
ious methods with 6-31G, cc-pVDZ, and cc-pVTZ basis sets. 
The dissociation energy D e was obtained with additional 
atomic calculations for the nitrogen atom. 



Re I^e D e 

A cm -1 kcal/mol 



6-31G 


PPT 


1 
1 




oO 


9908 


97 


1 «9i 
lOo. 


4:0 


iii 1 r 


n 

-u 




74 


4^9 


QA 


-1UO. 




CCSD 


-0. 


.006 


77 


77. 


53 


-6. 


.27 


CCSDT 


-0. 


.002 


06 


26. 


93 


-1 


.38 


CASSCF 


-0 


.003 


77 


29. 


35 


-8. 


.87 


CASPT2 


-0 


.000 


91 


3 


.71 


-3 


.85 


CASPT3 


-0. 


.001 


01 


8. 


.80 


-1. 


.62 


MR-CI 


-0 


.000 


33 


3 


.14 


-1 


.24 


MR-CI+Q 


0. 


.000 


10 


-0. 


.15 


-0. 


.02 


MR-ACPF 


-0 


.000 


16 


1. 


.68 


-0. 


.49 


MR-AQCC 


-0. 


.000 


20 


1. 


.99 


-0. 


.55 


L-CTSD(CU) 





.001 


11 


-15. 


.35 


-0. 


.67 


L-CTSD(MK) 





.000 


53 


-10 


.34 


-1. 


.21 


cc-pVDZ 


MR-CI+Q 


1. 


.120 


36 


2321. 


25 


200. 


59 


RHF 


-0 


.043 


06 


436. 


.76 


-88. 


.43 


CCSD 


-0. 


.007 


54 


87. 


.11 


-8. 


.62 


CCSDT 


-0. 


.001 


87 


24. 


90 


-1 


.79 


CASSCF 


-0. 


.005 


89 


43. 


95 


-3 


.58 


CASPT2 


-0. 


.001 


21 


4. 


.48 


-4. 


.11 


CASPT3 


-0 


.001 


07 


7. 


.77 


-2 


.40 


MR-CI 


-0 


.001 


10 


8. 


.71 


-2 


.83 


MR-ACPF 


-0 


.000 


47 


3. 


.34 


-0. 


.48 


MR-AQCC 


-0. 


.000 


62 


4. 


.39 


-0 


.48 


L-CTSD(CU) 


-0. 


.001 


25 


-4. 


.77 


-0. 


.42 


L-CTSD(MK) 


-0 


.002 


07 


-6. 


.20 


-1. 


.13 


cc-pVTZ 


MR-CI+Q 


1. 


.104 


76 


2332. 


40 


216. 


00 


RHF 


-0 


.037 


59 


398. 


68 


-95. 


61 


CCSD 


-0. 


.008 


05 


90. 


.91 


-8. 


.30 


CCSDT 


-0. 


.001 


65 


22. 


60 


-0 


.15 


CASSCF 


-0 


.001 


10 


18. 


63 


-12. 


24 


CASPT2 


-0. 


.000 


45 


-4. 


.75 


-6 


.99 


CASPT3 


-0 


.001 


23 


10. 


06 


-2 


.41 


MR-CI 


-0 


.001 


24 


10. 


58 


-4. 


.22 


MR-ACPF 


-0 


.000 


43 


3. 


.09 


-0. 


.37 


MR-AQCC 


-0. 


.000 


63 


4 


.87 


-0. 


.43 


L-CTSD(MK) 


-0 


.002 


49 


2. 


.02 


-1. 


.68 


exptl 


1. 


.107 


68 


2358. 


57 


228.4 



TABLE VIII: Size consistency errors {mE h ) of CISD, ACPF, 
AQCC, CCSD, L-CTSD(CU), and L-CTSD(MK) calcula- 
tions. 



Be + He Be + 2He Be + 3He Be + 4He 



CISD 


3.10 


6.23 


9.46 


12.83 


ACPF 


-0.56 


-0.76 


-0.86 


-0.92 


AQCC 


1.98 


2.38 


2.66 


2.91 


CCSD 


0.00 


0.00 


0.00 


0.00 


L-CTSD(CU) 


0.00 


0.00 


0.00 


0.00 


L-CTSD(MK) 


0.00 


0.00 


0.00 


0.00 




N 2 + He N 2 


+ 2He N 2 


+ 3He N 2 


+ 4He 


CISD 


1.96 


4.00 


6.12 


8.31 


ACPF 


-0.41 


-0.71 


-0.94 


-1.11 


AQCC 


-0.57 


-0.98 


-1.28 


-1.50 


CCSD 


0.00 


0.00 


0.00 


0.00 


L-CTSD(CU) 


0.00 


0.00 


0.00 


0.00 


L-CTSD(MK) 


0.00 


0.00 


0.00 


0.00 



ence methods, such as MR-ACPF and MR-AQCC, par- 
ticularly for energetic quantities such as the MAE and 
D e . The shapes of the curves in the intermediate regions 
looked somewhat like the CASPT3 curves, though with 
significantly smaller absolute errors. The spectroscopic 
constants uj e ,R e and the non-parallelity error from L- 
CTSD were slightly less accurate than from MR-ACPF 
and MR-AQCC and this was in part related to our nu- 
merical approximations in solving the CT equations. 

B. Size-consistency 

As is well recognized, size-consistency is a crucial re- 
quirement for any correlation model to obtain chemically 
accurate results in systems with many correlated elec- 
trons. As discussed in our initial work [4j, the L-CT 
theory is naturally size-consistent. One way to see this 
is to observe that the energy is obtained as the expec- 
tation value of an effective Hamiltonian that contains 
only connected contributions by virtue of its construc- 
tion via a commutator expansion Here we verify the 
size consistency property of L-CT theory through explicit 
numerical calculations on supermolecules. We have cho- 
sen to use supermolecules that contain more than one 
type of molecule since certain approximate size-extensive 
theories such as the ACPF and AQCC methods (which 
modify the non-size-consistent CISD method) are rigor- 
ously size-consistent only in the special case when the 
supermolecule is made of n noninteracting identical sub- 
systems. 

Table IVTlTl gives the size consistency errors of L-CTSD, 
CISD, CCSD, ACPF and AQCC calculations for the 
Be + n He and N2 + n He, respectively. All calcula- 
tions used the HF wavefunction as the reference and the 
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TABLE IX: Energy difference of all-electron and frozen-core atomic calculations, i.e. _E(all electron) — _E(frozen core), by various 
methods with 6-31G basis sets. 





Be 


No 


Ms 

o 


Ar 


Ca 


FCI 


-D ooosi 

yj • kjkjkjo J. 


-0 00078 
kj .kjkjkj t o 


-0 00976 


-0 001 Q5 

KJ • KJKJ -L UO 


-0 00399 

KJ • KJKJkjLiZj 


CISD 


-0 00075 
kj •\j\J\J i o 


-0 00077 
\j . kjkjkj t i 


-0 00959 

\J . KJKJ £i\J £i 


-0 001 87 

U • KJKJ -L O ( 


-0 009Q9 


\jVjui-/l i J 


-0 oooso 


-0 00078 


-0 00977 

KJ . KJKJ£i 1 I 


-0 001 Q6 

KJ .KJKJ -L UKJ 


-0 00399 

KJ .KJKJkj£i£i 


CCSD 


-0.00078 


-0.00077 


-0.00265 


-0.00189 


-0.00308 


AQCC 


-0.00156 


-0.00098 


-0.00452 


-0.00205 


-0.00505 


ACPF 


-0.00335 


-0.00090 


-0.00503 


-0.00199 


-0.00535 


MP3 


-0.00099 


-0.00076 


-0.00275 


-0.00185 


-0.00324 


MP2 


-0.00105 


-0.00091 


-0.00289 


-0.00219 


-0.00318 


L-CTSD(CU) 


0.05377 


0.11300 


0.03612 


0.03801 


0.02855 


L-CTSD(MK) 


0.00004 


0.00460 


-0.00289 


0.00660 


-0.00320 



molecules/atoms were each separated by a distance of 
1000 bohr. Size consistency implies the condition E(A + 
riB) = E(A) + nE(B). As can be seen, the L-CTSD and 
CCSD calculations are rigorously size-consistent while 
those of CISD, ACPF and AQCC steadily increase. 

Consider now the ACPF energy functional of the non- 
interacting system A + nB, given by 



r ACPF 
A+nB 



(H A )+n{H B ) 



1 + (2/N)(S A \5 A ) + (2n/N)(S B \S B ) 



(29) 



where (Ha) = (^aI-HaI^a), N is the total number of 
electrons, which is equal to Na + nNs, and 5 denotes 
the orthogonal correlation component of 'J, e.g. 'J'a = 
^oa + ^A- If A- = B in eqn. (|2"9"|) , we readily confirm that 
^(l+n)B = (^+ n )FB CPF an d the energy is size-consistent. 
The size consistency error in the functional is obtained 
as, 

e(n) 



pACPF 
- A+nB ' 



oACPF 
A 



n F, 



ACPF 



= RaRb(N a (H b ) + N B (H A )) - R 2 A N B (H B ) - R%N A (H A ) 
R A Rl+ R 2 A R B /n 

(30) 

where R A = N A + 2(S A \S A ) ■ The errors of ACPF and 
AQCC indeed appear to behave as the above function. 
(Note that e(n = oo) does not vanish). 



C. Density dependence 

Rather than considering the scaling behaviour of the 
energy as we increase the number of molecules, we can 
also consider the complementary trend of going to atoms 
with larger and larger nuclear charge Z (and conse- 
quently more and more electrons in the same region of 
space). In essence, this measures the density dependence 
of the energy. To study the behaviour of the CT and 
other methods under this condition, we chose five closed- 
shell atoms, two rare gas atoms (Ne[10e] and Ar[18e]) and 
three alkaline earth metals (Be[4e], Mg[12e] and Ca[20e]). 

Table llXl and figure [S] present the core-correlation ener- 
gies, defined as the energy difference between all-electron 



and frozen-core atomic calculations, using 6-31G basis 
sets. The frozen-core calculations correlate eight and two 
electrons in the valence orbitals for the rare gas atoms 
and alkaline earth metals, respectively, and represent 
98.3% (Be), 99.3% (Ne), 92.1% (Mg), 95.3% (Ar), and 
88.7% (Ca) of the all-electron correlation energies. The 
energy difference between the all-electron and frozen-core 
calculations is the core correlation energy from the core- 
valence and core-external excitations. Since we can re- 
gard valence-electron correlation as a size- intensive quan- 
tity that is described by a fixed, i.e. 0(1), small number 
of valence electrons at a roughly constant valence elec- 
tron density [57j , the rest of the correlation for the bulk 
of the electrons, i.e. the core correlation, must contain 
the main density dependence as we change the number 
of electrons and nuclear charge Z. 

Compared to the exact FCI core correlation energies, 
it is clear that ACPF, AQCC and L-CTSD(CU) have 
difficulty reproducing the correct behaviour. In par- 
ticular large errors are found in the ACPF and AQCC 
calculations for the alkaline earth metals and in the L- 
CTSD(CU) calculations of the rare gas atoms. By con- 
trast, the size-inconsistent CISD method as well as the 
MP2 and MP3 methods are able to capture the correct 
behaviour of the core correlation. This illustrates the dif- 
ficulty in finding an ad-hoc size-consistency correction, 
as employed in ACPF and AQCC, that works under all 
conditions. Most interestingly, the new operator decom- 
position in L-CTSD(MK) behaves much better than L- 
CTSD(CU) and reproduces the correct behaviour. 



D. FeO binding curve 

As a realistic example of a difficult multireference prob- 
lem, we calculated the potential curve for the ground 1 5 A 
state of the FeO molecule. ANO basis sets (U [H[ of 
DZP quality were used, [21sl5pl0d6/]/(5s4p3dl/) and 
[14s9p4<i]/(3s2plcf) for the Fe and O basis, respectively. 
To facilitate the setup for the multireference calculations, 
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TABLE X: Total energies of MR-CI+Q (E h ) and differences (mE b ) of various methods from MR-CI+Q for the ground 1 5 A 
state of FeO molecule. 





1.50A 


1.57A 


1.65 Jl 


1.72A 


2.00A 


MRPT-I-O 


1 V\7 (1^84^ 
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-LOO 1 .U 1 UU 1 


1 337 R730.9 
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73 
0. ( 


-4.12 


n ^7 

U.J/ 


O JT _L O 


41 7Q 


41 78 


4fl ^7 

*±U.O / 


38 4fi 


98 73 




91 78 


91 41 


on Q3 


on 37 


1 7 Q9 


MRACPF 


0.51 


0.40 


0.32 


0.29 


0.45 


MRAQCC 


4.61 


4.46 


4.30 


4.16 


3.79 


L-CTSD(CU) a 


c 


c 


9.17 


9.43 


7.67 


L-CTSD(CU) b 


c 


c 


6.85 


6.85 


3.51 


L-CTSD(MK) a 


3.09 


2.73 


3.00 


3.85 


3.52 


L-CTSD(MK) b 


0.83 


1.21 


1.66 


2.25 


-0.16 



a) t s = 3.0 x 10 _1 and r d = 5.0 x 10~ 2 . 

b) t s = 1.5 x 10" 1 and r d = 5.0 x 10" 2 . 

c) Not converged. 
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FIG. 6: Density scaling: Energy difference -E(all electron) 
^(frozen core) for atomic calculations shown in Tables ILXl 



the initial orbitals were obtained from closed-shell RHF 
calculations for the l 1 !] state. The ten lowest lying or- 
bitals for 20 electrons 

(lcT) 2 (2a) 2 (3cT) 2 (4 ( 7) 2 (5cT) 2 (6 ( 7) 2 (l7r) 4 (27r) 4 (31) 



were held frozen for the CASSCF and subsequent dy- 
namic correlation calculations. We verified that the er- 
rors made by this orbital restriction were almost constant 
within 1 raEh and thus would not affect the shapes of 
the potential curves. The orbital (7a) 2 was treated as an 
external core orbital, which was optimized by CASSCF 
and then correlated. The remaining 12 electrons were 
fully correlated with 12 active orbitals 

(8cT) 2 (9cT) 1 (10a) (llcT) (3^) 4 (4^) 2 (57r)°(l(5) 3 (32) 

(the occupations are based on the ROHF configuration of 
the 5 A state) for CAS denoted as (12e, 12o). This CAS 
is derived from Fe 3d and 4s orbitals, oxygen 2p orbitals, 
and the third bonding and antibonding ir orbitals, which 



TABLE XI: Spectroscopic constants for the ground 1 5 A state 
of FeO molecule. 



R e (A) u e (cm 1 ) 



CASSCF 


1.703 


691.1 


CASPT2 


1.620 


913.8 


CASPT3 


1.657 


755.2 


MR- CI 


1.641 


844.4 


MR-CI+Q 


1.635 


863.2 


MR-ACPF 


1.636 


863.5 


MR-AQCC 


1.637 


858.7 


L-CTSD(MK) a 


1.631 


914.0 


L-CTSD(MK) b 


1.630 


911.1 


exptl. 


1.616 


880 



a) 
b) 



T a = 3.0 x 10" 1 and r d = 5.0 x 
Ts = 1.5 x 10 _1 and T d = 5.0 x 



5.0 x 10" 2 . 
10~ 2 . 
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FIG. 7: Potential curve for the ground 1 5 A state of FeO molecule, a) r s = 3.0 x 10" 1 and r d = 5.0 x 10~ 2 . b) r s = 1.5 x 10 _1 
and r d = 5.0 x 10" 2 . 



are formed from the oxygen 2p a orbital mixing with some 
Fe 4p ff [Hil. 

Figure shows the potential curves of FeO computed 
by various multirefcrence methods. As exact energies 
are not available for this system, we report the differ- 
ences from MR-CI+Q energies in Table |Xj Clearly, 
the MR-ACPF, MR-AQCC, and L-CTSD(MK) curves 
are all very close to each other, while the MR-CI and 
CASPT2/CASPT3 curves are significantly further away. 
The MR-ACPF curve follows the MR-CI+Q curve with 
deviations of less than 0.5 m£"h, while the MR-AQCC 
curve is also nearly parallel with deviations of 3.8-4.6 
mE h . The L-CTSD(CU) and L-CTSD(MK) curves were 
shifted relative to each other; the L-CTSD(MK) ener- 
gies were significantly closer to the MRCI+Q energies, 
with deviations of less than 2.3 mEh. CASPT2 seemed 
to overestimate the correlation energy, while going to the 
third order CASPT3 over-corrected too much in the op- 
posite direction and strongly underestimated the corre- 
lation energy. 

Table IXII shows the spectroscopic constants measured 
from the potential curves. While the basis used is proba- 
bly too small for direct comparison to experiment, we see 
that in relation to the experimental results, MR-CI and 
related modifications MR-CI+Q, MR-ACPF and MR- 
AQCC give frequencies that are too low and bond-lengths 
that are too long, while L-CTSD(MK) gives frequencies 
that are too high and slightly improved bond-lengths. 



As we have already seen in the difference between the 
CASPT2 and CASPT3 curves, multireference perturba- 
tion theory seemed to break down in this molecule. 

E. Timings 

It is our intention that the CT theory should be practi- 
cally applicable to problems of reasonable size, and let us 
now examine the computational timings for the multiref- 
erence calculations on the FeO molecule we have just dis- 
cussed. These are shown in Tabic IXIII All timings were 
obtained on a single CPU of the Altex system (Itanium 
1.5GHz) at the Research Center for Computational Sci- 
ence, Okazaki. As can be seen, the MR-CI based meth- 
ods were two to three orders of magnitude more expen- 
sive than CASPT2. L-CTSD(MK) displayed very satis- 
factory performance. Even in our protoypte CT imple- 
mentation, which did not use point-group symmetry, the 
single-point energy calculation took less time than even 
the CASPT2 calculation, while providing a significantly 
better accuracy competitive with MR-ACPF. 

V. SUMMARY AND CONCLUSIONS 

We have been developing the Canonical Transforma- 
tion theory to describe dynamic correlation in multiref- 
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erence problems. The theory uses a size-extensive uni- 
tary exponential acting on a multireference function. In 
our initial work, we introduced a central approximation 
that rendered the manipulation of this ansatz practical, 
namely a cumulant-based operator decomposition. This 
choice of decomposition is not unique, however, and in 
the current work we introduced a new operator decompo- 
sition, based on the extended normal ordering of Mukher- 
jee and Kutzelnigg [2, [H, HH , which possesses attractive 
formal and conceptual features. 

We carried out calculations at the Linearised Canonon- 
ical Transformation Theory Singles and Doubles (L- 
CTSD) level using both our earlier cumulant-based and 
current Mukherjee-Kutzelnigg operator decompositions. 
In studies of the water, nitrogen, and iron-oxide binding 
curves, we found the accuracy of L-CTSD to be com- 
petitive with some of the best existing multireference 
methods such as the Multireference Averaged Coupled 
Pair Functional, while the computational cost was two- 
three orders of magnitude less and comparable to that of 
Complete- Active-Space Second Order Perturbation The- 
ory. Compared to our earlier work, our results and com- 
putational timings were greatly improved, in part due to 
the use of a new numerical algorithm for converging the 
Canonical Transformation equations. 



VII. APPENDIX: IMPLEMENTING 
CANONICAL TRANSFORMATION THEORY 

A. Recapitulation 

In our previous implementation of the CT algorithm 
[3| we solved the residual equations using the following 
skeletal algorithm 

1. Set up the electronic Hamiltonian H and the one- 
and two-particle density matrices of a reference 
wavefunction. 

2. Compute the transformed Hamiltonian (eqn. 



3. Compute the residuals of CT amplitude equations. 



iJ? = <[#i >3 ,a;-aJ]i j2 ) 
([ffi,2,a^-<Ji, 2 ) 



R p J 



(33) 
(34) 



4. Update the amplitudes by adding the precondi- 
tioned residuals 



A* <- AJ - B»/D> 

API <— API J?P4 I DPI 
A st ^~ A st ~ n st I u at 



(35) 
(36) 



where the factors 1 /D P and 1 /D p s 1 are the diagonal 
preconditioners. 
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TABLE XII: Timings for different multireference methods for 
a single point calculation on the FeO curve. Note that the 
L-CTSD calculation did not use point-group symmetry, while 
Civ symmetry was used in all the other calculations. The 
time for the CASSCF calculation is not included. 



Time/s 



CASPT2 

CASPT3 

MR-CI+Q 

MR-ACPF 

L-CTSD(MK) a 



5900 
17000 
158000 
168000 

4500 



a) The time for constructing density matrices is not 
included. 



In addition, we employed a somewhat complicated divi- 
sion of the optimisation process into different steps in- 
volving different classes of excitations in the A operator. 



B. Preconditioning and orthogonalisation 

Our primary concern in the current implementation 
was to improve the convergence of the CT equations. To 
achieve this, instead of using a diagonal preconditioner as 
in (|35[) . (|36p . we updated the amplitudes through an ex- 
act Newton step. The simplest way to define the Newton 
update is through the linear equation 



DP s; v v AAV y 

n pq, vw a avw 
^st.yz ^^yz 



-iff 



with 



DP, v 

s, y 

DPI, vw 

st, yz 



([[#1,2, 



a ?!]l,2, 



J 1,2) 



,P9 



]l,2) 
-st 1 



n 



1.2 



(37) 
(38) 



(39) 
(40) 



We can interpret the D matrices as the derivatives of 
the residual or Hessians of the energy. However, eqns. 
(|3"9")) . (|4TJ)l are non-optimal as the search directions (i.e. 
the components of A) within the first-order interacting 
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space, namely the singly-external, doubly-external, and 
semi-internal excitations 



(a ab - " ij 
(a ak 



~ <&)*o 



IJ 



(41) 
(42) 
(43) 



generate a non-orthogonal and even linearly-dependent 
basis. The large spread in eigenvalues of the overlap of 
the first-order interacting basis (|4*Tj) , (|4*2")) . (|4*3")) can then 
cause poor convergence of the linear equations (|37[) and 

To remedy this, we first orthogonalise the first-order 
interacting basis by diagonalising the overlap matrix S 
made of the one, two, and three-particle density matrices, 



(«-<)t( a «_ a j)) 



S. 



ijk, Iran 



7j* 



5' 



a6 



IJ71 

ilmk 



((«« -^(agf-oSi)) (a ^6) 
7« 



(44) 
(45) 
(46) 
(47) 



and change to the orthogonalised excitation operators a° 
and a" fc 



a o — 1/2/ a i \ . o~ 1/2/ afc \ 
a M = ~ °<J + u.ijk i a ij - a ak) 



(48) 

K b -4) (49) 



_ g-l/2 (n al, 



We can then solve the Newton equations (f3"T)l . (|3"8)) in 
this orthogonalised representation. To do so, the quan- 
tities A, R, and D are re-expressed in terms of a° and 



.,ab 



A = 




h A a J> af 


(50) 


K = 


([Sl,2, 


<-<]i, 2 ) 


(51) 


Rfi - 






(52) 


jja,b 

/■*) u 


([[Si* 


1,2, a«-<] 1)2 ) 


(53) 


f^ab, cd 

{A, V 


([[Sl,2 


, a i/ _ a cdJl,2, % - a abJl,2/ 


(54) 



The numbers of operators a" and a" are 0(a e) and 
0(a 2 e 2 ), respectively. Thus the additional cost of the 
transformation is 0(a 6 e) for the terms involving a° and 
0(a 4 e 2 ) for the terms involving a"'. The diagonalisation 
of the overlap matrix S for the semi-internal and singly- 
external (i.e. eqn. (j4~4"|) - (|l6"|) ) requires a cost of 0(a 9 ). 
While the scaling of these steps is relatively high, they are 
not expected to be a bottleneck for systems where con- 
ventional CASSCF calculations can be performed (for ex- 
ample, internally contracted CASPT2 also contains steps 



with such cost [3g, [37|])- However, if we were to use a 
large active space arising from e.g. a DMRG calculation, 
a different algorithm should be used. 

Let us consider now the condition number of D and 
the convergence characteristics of the Newton equations 
in the orthogonalised representation. If D were formed 
without any operator decomposition approximation, then 
D would represent the true Hessian of the energy with re- 
spect to an orthogonal set of directions in the first-order 
interacting space. The condition number of D would then 
be governed by the excitation energy between the refer- 
ence and excited states, which could be expected to be 
reasonable in most systems. Loosely speaking, we can 
regard the improved condition number of D as arising 
from the cancellation of small eigenvalues of D by the 
large eigenvalues of S^ 1 ^ 2 . However, such a cancella- 
tion is unstable, if we approximate D using the operator 
decomposition. Therefore to further improve the condi- 
tion number of D we discarded those operators a° and 
af which corresponded to small eigenvalues of S. The 
eigenvalue truncation thresholds are denoted hereafter as 
t s for the singly external and semi-internal and Td for the 
doubly external excitations. This limits the largest lin- 

l / 2 \ , n ,-i/*-\ 



ear combination amplitude coefficients (e.g. S J'') ap 



pearing in eqns. gSJ), Iggj) to 0(t s ' ) and 0{r d 1 ) 
respectively, preserving numerical stability in the am- 
plitude equations. Typically, we used r s < 1CP 1 and 
Td < 10~ 2 . These cutoffs appear large because of the 
extreme degeneracy of the first-order interacting space 
near equilibrium, and because of the incomplete removal 
of the poorly conditioned components, due to the slight 
incompatibility (unstable cancellation) between the ap- 
proximate Hessian and the overlap matrix in this space. 
Linear dependency is particularly strong near equilibrium 
because some of the active orbitals which are being ex- 
cited by A have nearly zero occupancy. However, the 
contribution of the neglected terms to the energy is small. 
For N2 with the cc-pVDZ basis, the size of the effective 
orthogonalised operator space was 410 (Rnn — 1-0), 698 
(Rnn — 1.6), and 938 (Rnn = 3.0) indicating that over 
50% of the operator basis was truncated near equilib- 
rium. In the dissociation region of the potential energy 
curves studied here, truncation did not occur. 

In our previous work, we encountered numerical dif- 
ficulties in using singly external excitation operators in 
conjunction with doubles, i.e. for L-CTSD. We now see 
that the reason is the linear dependency between singly 
external and semi-internal excitations, which appears as 
non-zero overlap in S^jki (eqn. (|45|) ). The orthogonal- 
isation fixes this issue, and thus we have used L-CTSD 
as the standard L-CT model in this work. This should 
be naturally superior to L-CTD as it includes orbital re- 
laxation and extra correlation such as three- or higher- 
particle excitations from the direct product of singles and 
doubles. 

Using the Newton update as described above, we ob- 
served efficient convergence in the CT amplitude equa- 
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TABLE XIII: The total energies (E b ) of L-CTSD(MK) for 
the bond breaking of the N2 molecule with CAS(6e,6o) and 
various basis sets using the exact and approximate (cumulant) 
three-particle density matrices for orthogonalisation. 





lA 


2A 


3l 


6-31G 








exact orthog. 


-109.04585 


-108.86156 


-108.84155 


cumulant orthog 


-109.04566 


-108.86220 


-108.84155 


diff (mEh) 


+0.19 


-0.64 


0.0 


cc-pVDZ 








exact orthog. 


-109.22774 


-108.98214 


-108.96009 


cumulant orthog 


-109.22752 


-108.98305 


-108.96009 


diff (mE h ) 


+0.22 


-0.91 


0.0 


cc-pVTZ 








exact orthog. 


-109.33525 


-109.05709 


-109.03011 


cumulant orthog 


not conv. 


-109.05824 


-109.03011 


diff (mE h ) 




-1.15 


0.0 



tions. Typically only 10 Newton steps would be required 
to converge the amplitudes in multireference calcula- 
tions. Convergence behavior of the amplitude equations 
in L-CTSD(CU) and L-CTSD(MK) was generally simi- 
lar, but there were some cases where we could converge 
the L-CTSD(MK) but not the L-CTSD(CU) calculations 



with the standard truncation thresholds, for example at 
R NN = l& for N 2 (cc-pVTZ), as discussed in Sec. IIVAI 



C. Operator orthogonalisation with cumulant 
density matrix 

Rather than using the exact three-particle density ma- 
trix for the orthogonalisation procedure described above, 
we could also imagine using its cumulant decomposition 
in terms of the one- and two-particle density matrices, 
eqn. (flU)l . 

Tabic IXIIII shows the differences of the total ener- 
gies computed using the amplitude operators that are 
orthogonalised with the exact and approximate (cumu- 
lant) three-particle density matrices for the N2 potential 
energy curves discussed in section ITV Al With the trun- 
cation threshold r s = 1CT 1 and = 10 -2 and using var- 
ious basis sets, the energies from both orthogonalisations 
were generally in good agreement within a deviation of 
1.2 mi?h- However, at Rnn — 1.0A with the cc-pVTZ 
basis set, the L-CTSD(MK) calculation with the cumu- 
lant based orthogonalisation did not converge. 
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